BeginPackage["Backpropagation`"]

sigmoid::usage = "sigmoid[x,opts___Rule]"
bpnTest::usage = "bpnTest[hiddenWts,outputWts,ioPairVectors,opts___Rule]"
bpnStandard::usage = "bpnStandard[inNumber, hidNumber, outNumber,ioPairs, eta, numIters]"
bpnBias::usage = "bpnBias[inNumber, hidNumber, outNumber,ioPairs, eta, numIters]"
bpnMomentum::usage = "bpnMomentum[inNumber,hidNumber,outNumber,ioPairs,eta,
     alpha,numIters]"
bpnMomentumSmart::usage = "bpnMomentumSmart[inNumber,hidNumber,outNumber,ioPairs,eta,
				alpha,numIters]"
bpnCompete::usage = "bpnCompete[inNumber,hidNumber,outNumber,ioPairs,eta,numIters]"


Begin["`Private`"]    (* begin the private context *)


Options[sigmoid] = {xShift->0,yShift->0,temperature->1};
Options[bpnTest] = {printAll->False,bias->False};
 
sigmoid[x_,opts___Rule] :=
	Module[{xshft,yshft,temp},
		xshft = xShift /. {opts} /. Options[sigmoid];
		yshft = yShift /. {opts} /. Options[sigmoid];
		temp = temperature /. {opts} /. Options[sigmoid];
		yshft+1/(1+E^(-(x-xshft)/temp)) //N
		]
 
bpnTest[hiddenWts_,outputWts_,ioPairVectors_,opts___Rule] :=
  Module[{inputs,hidden,outputs,desired,errors,i,len=Length[inputs],
     prntAll,errorTotal,errorSum,bias},
   prntAll= printAll /. {opts} /. Options[bpnTest];
   biasVal = bias /. {opts} /. Options[bpnTest];
   inputs=Map[First,ioPairVectors];
   If[biasVal,inputs=Map[Append[#,1.0]&,inputs] ];
   desired=Map[Last,ioPairVectors];
   hidden=sigmoid[inputs.Transpose[hiddenWts]];
   If[biasVal,hidden = Map[Append[#,1.0]&,hidden] ];
   outputs=sigmoid[hidden.Transpose[outputWts]];
   errors= desired-outputs;
   If[prntAll,Print["ioPairs:"];Print[ ];Print[ioPairVectors];
   				Print[ ];Print["inputs:"];Print[ ];Print[inputs];
   				Print[ ];Print["hidden-layer outputs:"];
   				Print[hidden];Print[ ];
   				Print["output-layer outputs:"];Print[ ];
   				Print[outputs];Print[ ];Print["errors:"];
   				Print[errors];Print[ ];  ]; (* end of If *)
   For[i=1,i<=len,i++,Print[" Output ",i," = ",outputs[[i]]," desired = ",
		   desired[[i]]," Error = ",errors[[i]]];Print[];  ];			(* end of For *)
   errorSum = Apply[Plus,errors^2,2]; (* second level *)
   errorTotal = Apply[Plus,errorSum];
   Print["Mean Squared Error = ",errorTotal/len];
   ]					(* end of Module *)
 

bpnStandard[inNumber_, hidNumber_, outNumber_,ioPairs_, eta_, numIters_] :=
  Module[{errors,hidWts,outWts,ioP,inputs,outDesired,hidOuts,
  			outputs, outErrors,outDelta,hidDelta},
	hidWts = Table[Table[Random[Real,{-0.1,0.1}],{inNumber}],{hidNumber}];
	outWts = Table[Table[Random[Real,{-0.1,0.1}],{hidNumber}],{outNumber}];
	errors = Table[
			(* select ioPair *)										
	  ioP=ioPairs[[Random[Integer,{1,Length[ioPairs]}]]];
	  inputs=ioP[[1]];
	  outDesired=ioP[[2]];
	  		(* forward pass *)
	  hidOuts = sigmoid[hidWts.inputs];
	  outputs = sigmoid[outWts.hidOuts];
	  		(* determine errors and deltas *)
	  outErrors = outDesired-outputs;
	  outDelta= outErrors (outputs (1-outputs));
	  hidDelta=(hidOuts (1-hidOuts)) Transpose[outWts].outDelta;
	  		(* update weights *)
	  outWts += eta Outer[Times,outDelta,hidOuts];
	  hidWts += eta Outer[Times,hidDelta,inputs];
	  		(* add squared error to Table *)
	  outErrors.outErrors,{numIters}];  (* end of Table *)
	Return[{hidWts,outWts,errors}];
	];									(* end of Module *)

 

bpnBias[inNumber_, hidNumber_, outNumber_,ioPairs_, eta_, numIters_] :=
  Module[{errors,hidWts,outWts,ioP,inputs,outDesired,hidOuts,
  			outputs, outErrors,outDelta,hidDelta},
	hidWts = Table[Table[Random[Real,{-0.1,0.1}],{inNumber+1}],{hidNumber}];
	outWts = Table[Table[Random[Real,{-0.1,0.1}],{hidNumber+1}], {outNumber}];
	errorList = Table[
			(* select ioPair *)										
	  ioP=ioPairs[[Random[Integer,{1,Length[ioPairs]}]]];
	  inputs=Append[ioP[[1]],1.0]; (* bias mod *)
	  outDesired=ioP[[2]];
	  		(* forward pass *)
	  hidOuts = sigmoid[hidWts.inputs];
	  outInputs = Append[hidOuts,1.0];  (* bias mod *)
	  outputs = sigmoid[outWts.outInputs];
	  		(* determine errors and deltas *)
	  outErrors = outDesired-outputs;
	  outDelta= outErrors (outputs (1-outputs));
	  hidDelta=(outInputs (1-outInputs)) * Transpose[outWts].outDelta;
	  		(* update weights *)
	  outWts += eta Outer[Times,outDelta,outInputs];
	  hidWts += eta Drop[Outer[Times,hidDelta,inputs],-1];  (* bias mod *)
	  		(* add squared error to Table *)
	  outErrors.outErrors,{numIters}];  (* end of Table *)
	  Print["New hidden-layer weight matrix: "];
  Print[]; Print[hidWts];Print[];
  Print["New output-layer weight matrix: "];
  Print[]; Print[outWts];Print[];  
  bpnTest[hidWts,outWts,ioPairs];   (* check how close we are *)
  errorPlot = ListPlot[errorList, PlotJoined->True];
  Return[{hidWts,outWts,errorList,errorPlot}];
	];									(* end of Module *)
 

bpnMomentum[inNumber_,hidNumber_,outNumber_,ioPairs_,eta_,
     alpha_,numIters_] :=
  Module[{hidWts,outWts,ioP,inputs,hidOuts,outputs,outDesired, 
    hidLastDelta,outLastDelta,outDelta,hidDelta,outErrors},
    hidWts = Table[Table[Random[Real,{-0.5,0.5}],{inNumber}],{hidNumber}];
    outWts = Table[Table[Random[Real,{-0.5,0.5}],{hidNumber}],{outNumber}];
    hidLastDelta = Table[Table[0,{inNumber}],{hidNumber}];
    outLastDelta = Table[Table[0,{hidNumber}],{outNumber}];
    errorList = Table[
					         (* begin forward pass *)
	     ioP=ioPairs[[Random[Integer,{1,Length[ioPairs]}]]];
	     inputs=ioP[[1]];
     	outDesired=ioP[[2]];
	     hidOuts = sigmoid[hidWts.inputs];  (* hidden-layer outputs *)
	     outputs = sigmoid[outWts.hidOuts]; (* output-layer outputs *)
					        (* calculate errors *)
	     outErrors = outDesired-outputs;
	     outDelta= outErrors (outputs (1-outputs));
	     hidDelta=(hidOuts (1-hidOuts)) Transpose[outWts].outDelta;
					      (* update weights *)
	     outLastDelta= eta Outer[Times,outDelta,hidOuts]+alpha outLastDelta;
	     outWts += outLastDelta;
	     hidLastDelta = eta Outer[Times,hidDelta,inputs]+
						   	alpha hidLastDelta;
	     hidWts += hidLastDelta;
		    outErrors.outErrors, (* this puts the error on the list *)
		    {numIters}]	;	 (* this many times, Table ends here *)
  Print["New hidden-layer weight matrix: "];
  Print[]; Print[hidWts];Print[];
  Print["New output-layer weight matrix: "];
  Print[]; Print[outWts];Print[];  
  bpnTest[hidWts,outWts,ioPairs,bias->False,printAll->False];
  errorPlot = ListPlot[errorList, PlotJoined->True];
  Return[{hidWts,outWts,errorList,errorPlot}];
  ]						(* end of Module *)
 

bpnMomentumSmart[inNumber_,hidNumber_,outNumber_,ioPairs_,eta_,
				alpha_,numIters_] :=
  Module[{hidWts,outWts,ioP,inputs,hidOuts,outputs,outDesired, 
			 hidLastDelta,outLastDelta,outDelta,hidDelta,outErrors},
    hidWts = Table[Table[Random[Real,{-0.5,0.5}],{inNumber}],{hidNumber}];
    outWts = Table[Table[Random[Real,{-0.5,0.5}],{hidNumber}],{outNumber}];
    hidLastDelta = Table[Table[0,{inNumber}],{hidNumber}];
    outLastDelta = Table[Table[0,{hidNumber}],{outNumber}];
    errorList = Table[
					      (* begin forward pass *)
	    ioP=ioPairs[[Random[Integer,{1,Length[ioPairs]}]]];
	    inputs=ioP[[1]];
	    outDesired=ioP[[2]];
	    hidOuts = sigmoid[hidWts.inputs];  (* hidden-layer outputs *)
    	outputs = sigmoid[outWts.hidOuts]; (* output-layer outputs *)
				    	(* calculate errors *)
	    outErrors = outDesired-outputs;
	    If[First[Abs[outErrors]]>0.1,
	      outDelta= outErrors (outputs (1-outputs));
	      hidDelta=(hidOuts (1-hidOuts)) Transpose[outWts].outDelta;
					      (* update weights *)
	      outLastDelta= eta Outer[Times,outDelta,hidOuts]+
							    alpha outLastDelta;
      	outWts += outLastDelta;
	      hidLastDelta = eta Outer[Times,hidDelta,inputs]+
						       	alpha hidLastDelta;
	      hidWts += hidLastDelta,Continue]; (* end of If *)
		  outErrors.outErrors, (* this puts the error on the list *)
		  {numIters}]	;	 (* this many times, Table ends here *)
  Print["New hidden-layer weight matrix: "];
  Print[]; Print[hidWts];Print[];
  Print["New output-layer weight matrix: "];
  Print[]; Print[outWts];Print[];  
  bpnTest[hidWts,outWts,ioPairs,bias->False,printAll->False];
  errorPlot = ListPlot[errorList, PlotJoined->True];
  Return[{hidWts,outWts,errorList,errorPlot}];
  ]						(* end of Module *)
 

bpnCompete[inNumber_,hidNumber_,outNumber_,ioPairs_,eta_,numIters_] :=
 Module[{hidWts,outWts,ioP,inputs,hidOuts,outputs,outDesired, 
			 outInputs,hidEps,outEps,outDelta,hidPos, outPos, hidDelta,outErrors},
  hidWts = Table[Table[Random[Real,{-0.5,0.5}],{inNumber}]{hidNumber}];
  outWts = Table[Table[Random[Real,{-0.5,0.5}],{hidNumber}],{outNumber}];
  errorList = Table[    (* begin forward pass *)
  ioP=ioPairs[[Random[Integer,{1,Length[ioPairs]}]]];
	 inputs=ioP[[1]]; 
	 outDesired=ioP[[2]];
	 hidOuts = sigmoid[hidWts.inputs];
	 outputs = sigmoid[outWts.hidOuts]; 
	 outErrors = outDesired-outputs;      (* calculate errors *)
	 outDelta= outErrors (outputs (1-outputs));
 	hidDelta=(hidOuts (1-hidOuts)) Transpose[outWts].outDelta;
					    (* index of max delta *)
 	outPos = First[Flatten[Position[Abs[outDelta],Max[Abs[outDelta]]]]];
	 outEps = outDelta[[outPos]];  (* max value *)
	 outDelta=Table[-1/4 outEps,{Length[outDelta]}]; (* new outDelta table *)
	 outDelta[[outPos]] = outEps;  (* reset this one  *)
				    (* index of max delta *)
	 hidPos = First[Flatten[Position[Abs[hidDelta],Max[Abs[hidDelta]]]]];
	 hidEps = hidDelta[[hidPos]];  (* max value *)
	 hidDelta=Table[-1/4 hidEps,{Length[hidDelta]}]; (* new outDelta table *)
		hidDelta[[hidPos]] = hidEps;  (* reset this one  *)
	 outWts +=eta Outer[Times,outDelta,hidOuts];
	 hidWts += eta Outer[Times,hidDelta,inputs];
		outErrors.outErrors, (* this puts the error on the list *)
		{numIters}]	;	 (* this many times, Table ends here *)
 Print["New hidden-layer weight matrix: "];
 Print[ ]; Print[hidWts];Print[ ];
 Print["New output-layer weight matrix: "];
 Print[ ]; Print[outWts];Print[ ];  
 bpnTest[hidWts,outWts,ioPairs,bias->False,printAll->False]; 
 errorPlot = ListPlot[errorList, PlotJoined->True];
 Return[{hidWts,outWts,errorList,errorPlot}];
 ]						(* end of Module *)	
 


End[]         (* end the private context *)

EndPackage[]  (* end the package context *)